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L  INTRODUCTION. 


Over  the  past  several  years  extensive  efforts  have  been  made  in  using  adaptive 
strategies  to  solve  partial  differential  equations  [2,  3]*.  In  this  report,  we  consider  a 
local  mesh  refinement  procedure  for  two^limensional  parabolic  partial  differential  sys¬ 
tems  where  fine  meshes  are  introduced  in  regions  where  greater  resolution  is  deemed 
necessary.  Our  approach  permits  finer  meshes  to  overlap  elements  of  coarser  ones 
and  is  related  to  an  earlier  effort  on  h-refinement  methods  for  one-dimensional  para¬ 
bolic  problems  [5,  7,  10]. 

We  consider  an  initial-boundary  value  problem  for  an  m-dimensional  vector  sys¬ 
tem  having  the  form 

u,  +  f(x  O’  .u,u,  .Uy )  =  [Di(x  ,y  ,f  .u)Ujt ],  +  [Djfx  o’  ,t  ,u)Uv  ]y,  (x  ,y )  e  O,  t  >  0, 


(la) 

u(xo’,0)  =  Uo(xo’), 

(x,y)6  0vj30, 

(lb) 

u(xo',0  =  go(*.y.O. 

(x,y)GdQp,  r>0. 

(Ic) 

-»■  D2UyV  =  g/v0fO’^)>  (x0')€9Ov,  t  >  0. 

(Id) 

The  domain  O  is  the  rectangle  {(xo’) 

la<x<b,c<y<d}  with 

boundary 

30  =  BQq  and  unit  outer  normal  T|  :=  [t|^T|^]^.  The  system  in  Eq.  (1)  is 

assumed  to  be  well  posed  and  parabolic,  i.e.,  Dj  and  D2  are  positive  definite.  We  do 
not  expect  that  our  methods  will  be  able  to  solve  all  problems  having  this  generality, 
but  our  one-dimensional  procedure  [10]  has  worked  well  on  a  wide  range  of  linear  and 
nonlinear  problems. 


*  References  ere  listed  is  the  end  of  this  repoiL 


Our  approach  begins  with  the  solution  of  Eq.  (1)  on  a  uniform  space-time  grid 
using  finite  elements  in  space  and  the  backward  Euler  method  in  time.  At  the  end  of 
each  time  step,  an  indication  of  the  local  discretization  error  is  generated  on  each  finite 
element  In  our  initial  investigation  of  one-dimensional  problems  [5,  7],  we  used  an 
h-refinement  (Richardson’s  extrapolation)  procedure  to  compute  a  local  error  indicator. 
This  has  subsequently  been  abandoned  in  favor  of  a  p-refinement  approach  [10],  which 
increases  the  order  of  the  trial  space  instead  of  reducing  the  mesh  spacing.  The  p- 
refinement  strategy  employs  nodal  superconvergence  to  improve  computational 
efficiency  and  it  can  be  used  to  generate  an  asymptotically  correct  estimate  of  the 
discretization  error  [1,  10].  Elements  having  high  error  are  grouped  into  rectangular 
regions  called  megagrids  using  a  nearest  neighbor  clustering  algorithm  (cf.  Berger  and 
Oliger  [4]).  Overlapping  fine  uniform  grids  are  generated  within  the  megagrids  and 
Eq.  (1)  is  solved  again  on  these  grids.  This  process  is  repeated  until  a  prescribed  local 
error  tolerance  is  satisfied.  An  illustration  of  a  coarse  spatial  mesh  with  two 
megagrids  and  three  fine  grids  is  shown  in  Figure  1. 

A  tree  is  a  natural  data  structure  to  manage  the  information  associated  with  all  of 
the  grids.  Nodes  of  the  tree  represent  data  at  the  megagrid  level,  with  finer  megagrids 
regarded  as  offspring  of  coarser  ones.  Information  associated  with  overlapping  fine 
grids  within  each  megagrid  is  stored  as  records  at  the  nodes  of  the  tree. 

A  finite  element  problem  is  formulated  and  solved  on  each  grid  within  a 
megagrid.  This  necessitates  the  prescription  of  appropriate  initial  and  boundary  condi¬ 
tions  on  each  space-time  grid.  Since  our  temporal  integration  is  implicit,  prescribing 
boundary  conditions  is  particularly  complex  in  regions  where  meshes  overlap  (Figure 
1).  An  iterative  procedure,  analogous  to  Schwarz  alternation  (cf.  Dihn  et  al.  [6]),  is 
used  to  successively  calculate  solutions  on  fine  grids  within  each  megagrid.  We 


Figure  1.  Coarse  spadal  background  mesh  with  two  offspring  megagrids 
(marked  with  diamonds  and  squares)  and  their  local  grids.  High-eiror  ele¬ 
ments  of  the  coarse  mesh  are  indicated  by  x’s. 
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observe  that  this  procedure  converges  for  a  variety  of  problems,  but  have  no  analysis 
demonstrating  either  convergence  or  stability.  Starius  [11]  obtained  some  stability 
results  on  a  similar  method  for  hyperbolic  equations. 

A  description  of  the  data  structures  and  the  local  refinement  procedure  is  given  in 
Section  n.  In  Section  m  we  present  the  finite  element  method  and  the  local  error  esti¬ 
mation  technique.  Section  IV  contains  some  prelinunary  computation  results  on  three 
linear  parabolic  problems.  Our  conclusions  and  plans  for  further  improvements  are 
described  in  Section  V.  The  examples  indicate  that  the  error  estimation  procedure 
converges  to  the  true  discretization  error  as  the  mesh  is  refined,  and  the  solution  pro¬ 
cedure  based  on  the  Schwarz  alternating  technique  converges. 

n.  LOCAL  REFINEMENT  AND  DATA  STRUCTURES. 

We  outline  our  procedure  for  solving  Eq.  (1)  on  an  arbitrary  hexahedral  megagrid 
The  domain  (O  :=  {(x,y)  I  a<x  <  P,  y<y  <  5);  p  and  are  the 
times  at  the  beginning  and  end  of  the  time  step,  respectively;  F  and  S  point  to  the 
parent  and  offspring  megagrids,  respectively;  and  L  is  the  record  of  information  for 
the  a  local  rectangular  grids  within  R . 

A  top  level  description  of  our  local  refinement  algorithm  is  presented  in  Figure  2. 
Solution  and  error  indicators  are  generated  on  Jt  using  procedure  solve  (Figure  3). 
Elements  where  the  error  indicator  exceeds  a  prescribed  tolerance  to/  are  partitioned 
into  rectangular  regions  using  the  nearest  neighbor  clustering  algorithm.  As  noted,  we 
call  these  regions  megagrids.  Berger  and  Oliger’s  [4]  bisection  and  merging  procedure 
is  used  to  generate  local  uniform  fine  grids  for  each  megagrid.  Local  grids  within  a 
megagrid  can  overlap,  but  the  megagrids  are  independent  of  each  other,  hence,  each 
offspring  megagrid  may  have  different  spatial  and  temporal  refinement  factors.  This 


also  reduces  communication  between  megagrids,  and  thus,  simplifies  the  computation 
of  initial  conditions  on  offspring  megagrids.  This  representation  may  additionally  be 
suitable  for  execution  on  parallel  computers.  Temporal  refinement  factors  are  calcu¬ 
lated  and  solutions  are  recursively  generated  for  each  megagrid. 


procedure  locref  (R  (to^  ,<7  ,5  X  ),to/ ); 

begin 
solve 

if  any  error  indicator  >  tol  then 
be^n 

Form  offspring  megagrids; 
for  j  :=  1  to  number  of  offspring  do 
Qeate  local  rectangular  grids; 
for  j  :=  1  to  number  of  offspring  do 
Calculate  the  temporal  refeement  factor  tref  [j]; 
for  j  ;=  1  to  numbCT  of  offspring  do 
for  I  :*  1  to  tref  [y  ]  do 
begin 

p[i]  ;=  p+(i-l)*(q-p)/treflj]; 
qti]  :=  p(i]  +(q-p)/trcf[j]; 
locref  (R  (©[y  [1  ],q  [1  ], 

R  (©^  X  ,5  X  )X  [y  ]X  [y  ]),ro/ ) 

end 

end 

end  {  locref  }; 


Figure  2.  Recursive  local  refinement  algorithm  for  the  solution  of  Eq.  (1)  on 
R  (fuqj  ,<7  X  X  X )  with  an  error  tolerance  tol . 


In  order  to  solve  the  problem  in  Eq.  (1),  the  procedure  locref  is  invoked  on  the 
coarse  grids  k  =  0,  1,  •  •  • .  Solutions  satisfying  the  prescribed 

accuracy  requirements  are  generated  at  each  time  tt,k  =  \,1,  ■  -  ■  . 

The  solution  on  a  megagrid  Ri(ii,p  ,qJF ^  JL)  is  described  by  the  procedure  solve 
of  Figure  3.  Initial  conditions  are  generated  for  each  local  computation  grid  contained 
in  R.  Following  this,  we  compute  an  initial  guess  for  the  boundary  conditions  of  the 


local  grids  using  either  the  prescribed  boundary  data  at  physical  boundaries  or  linear 
interpolation  in  time  from  the  parent  megagrid  of  R.  A  finite  element  solution  is  gen¬ 
erated  for  one  of  the  local  grids  and  its  solution  is  used  to  update  boundary  conditions 
on  all  other  intersecting  local  grids.  This  solution  process  is  repeated  on  each  local 
grid  in  turn  until  satisfactory  convergence  is  attained.  Our  procedure  is,  thus,  similar 
to  the  Schwarz  alternating  principle  for  elliptic  problems,  which  has  been  used  recently 
to  develop  domain  decomposition  methods  for  parallel  computation  [6,  8]. 


procedure  solve  {R (o)^ 
begin 

for  i  :s  1  to  number  of  local  grids  do 
begin 

Compute  initial  conditions  for  local  grid 

^  (C^m  )«'  »(ym  )» )i  )»■  )» 

Compute  boundary  conditions  for 
T  (ix„  )i  ,(y^)i  ,(d,  )i  Xdy  )i  ); 

end 

for  j  :s  1  to  number  of  iterations  do 
for  i  :=  1  to  number  of  local  grids  do 
begin 

Solve  the  finite  element  problem  for  Eq.  (1)  on 
T  ((^m  )i  )i  .(d,  )i  ,(^  )i  ); 

if  j  =  number  of  iterations  then 
Compute  error  on  T  ),•  ,(y^  ),•  ,(dx  )i  .(dy  ),•  ,J,- ) 
Update  appropriate  boundary  conditions 
end 

end  {  solve  }  ; 


Figure  3.  Solution  algorithm  on  megagrid  R {(0,p . 


A  local  grid  is  denoted  as  T(Xn,y^4x4yyS%  Each  local  rectangular  grid  is 
characterized  by  the  coordinates  of  its  center  the  lengths  of  its  sides  and 

dy ,  and  the  slope  s  of  a  side  of  the  rectangle.  In  order  to  avoid  ambiguity,  we  choose 
J  ^  0  and  let  d,  correspond  to  this  side  (Figure  1).  The  number  of  elements  mi  xni 


on  local  grid  T,-  =  T((Xn)i,(yn)i,(dx)i,(dy)iA)  **  determined  by  a  single  mesh  spacing 
parameter  hn  as  m,-  =  round (dj^/hn)  and  /ij  -  round (d^lhn).  Thus,  each  local  grid  in 
R  has  approximately  the  same  spadal  resolution.  Many  details  of  this  algorithm  have 
been  omitted  and  addidonal  informadon  is  presented  in  Moore  [9].  For  example,  a 
strategy  has  been  developed  for  storing  the  finite  element  soludon  at  p  and  q  without 
unnecessary  duplicadon  or  copying  of  informadon. 

Inidal  condidons  for  each  local  grid  are  either  determined  from  Eq.  (lb)  when 
p  =  0  or  by  bilinear  interpoladon  using  the  finest  grids  available  in  the  tree  structure  at 
time  p  >  0.  Isolating  local  grids  within  megagrids  gready  simplifies  the  search  for 
data  needed  for  this  bilinear  interpoladon.  Thus,  the  search  for  a  soludon  value  at  an 
arbitrary  point  is  performed  at  the  megagrid  level  until  the  finest  megagrid  containing 
the  point  has  been  identified.  The  local  grids  of  this  finest  megagrid  provide  the 
necessary  interpolation  data.  Scanning  the  points  of  a  grid  in  a  predetermined  order 
can  be  used  to  further  reduce  the  complexity  of  the  search  procedure. 

Similar  considerations  are  required  to  determine  boundary  conditions  on  grid 
edges  that  are  not  subsets  of  3£2.  Our  one-dimensional  techniques  [10]  and  the  expli¬ 
cit  finite  difference  procedures  of  Berger  and  Oliger  [4]  used  the  notion  of  a  "buffer" 
to  apply  boundary  conditions.  The  idea  is  to  enlarge  a  local  rectangular  grid  by 
increasing  and  dy  by  two  or  four  elements  so  that  "artificial  boundary  conditions" 
may  be  obtained  from  data  in  low-error  regions.  However,  in  regions  where  local 
grids  overlap,  accurate  boundary  conditions  cannot  be  obtained  from  parent  grid  data 
even  with  a  buffer.  Buffers  do  provide  accurate  boundary  conditions  in  regions  where 
grids  do  not  overlap  and,  for  this  reason,  we  continue  to  use  them. 


8 

Dirichlet  boundary  conditions  are  obtained  on  the  edges  of  buffered  local  grids  by 
piecewise  bilinear  interpolation  in  time  using  solution  values  from  the  parent  megagrid. 
In  non-overlapping  buffered  regions,  the  interpolated  boundary  conditions  satisfy  the 
prescribed  error  tolerance,  and  are  thus  expected  to  produce  acceptable  accuracy.  As 
noted,  accurate  boundary  conditions  are  obtained  in  regions  where  local  grids  overlap 
by  means  of  the  Schwarz  alternating  principle.  Hence,  we  iiutially  solve  a  finite  ele¬ 
ment  problem  on  local  grid  7  ^  of  R,  realizing  that  the  interpolated  boundary  data  may 
be  inaccurate  in  regions  where  T  i  intersects  other  local  grids.  In  solving  the  problem 
on  T2  we  use  boundary  data  from  T  j  with  bilinear  interpolation  in  regions  where  T  j 
and  T 2  intersect  This  sequential  updating  procedure  can  be  continued  iteratively  until 
satisfactory  convergence  is  obtained.  In  practice,  we  halt  the  iteration  after  a  few 
cycles  and  compute  an  error  estimate  for  each  local  grid  in  R.  The  grids  of  R  are 
refined  if  the  error  tolerance  is  not  satisfied.  Thus,  we  do  not  distinguish  between 
failure  of  the  Schwarz  iteration  to  converge  and  failure  to  satisfy  prescribed  accuracy 
conditions. 

Treatment  of  situations  'where  local  grids  overlap  are  considerably  more  com¬ 
plex.  A  second  complication  arises  when  a  local  grid  crosses  the  boundary  of 
o, 

where  the  subscript  F  denotes  the  parent  megagrid  of  R.  These  issues  are 
i=l 

handled  by  regridding  as  described  in  Moore  [9]. 

in.  SPATUL  AND  TEMPORAL  DISCRETIZATION, 

As  noted,  the  partial  differential  system  in  Eq.  (1)  is  discretized  on  a  local  grid  T 
of  R  using  a  finite  element  Galerkin  procedure  in  space  and  the  backward  Euler 
method  in  time.  For  each  time  /  e  [p,^],  we  assume  that  ue  Hg  and  select  a  test 


function  ve  //q  ,  where  denotes  the  usual  Sobolev  space.  Functions  that  further 
satisfy  Diiichlet  conditions  on  97  are  said  to  belong  to  H£,  while  functions  satisfying 
trivial  Dirichlet  conditions  belong  io  Hq. 

The  Galeikin  form  of  Eq.  (1)  on  T  is 


(v,u, )  +  (v,f(-,-^  ,u.u,  ,Uy ))  +  i4  (v.u)  =  j  ^iffds ,  for  all  V  €  //  J  .  (2a) 


where 


(v.u)  =  I  udxdy ,  (2b) 

A  (v,u)  =  I  [vjDiCx  ,y  a  ,u)u,  +  Vy^D2(x  ^  .f  ,u)iiy  ]dxdy.  (2c) 

Initial  conditions  are  required  at  p  =  0  and  these  can  be  obtained,  e.g.,  by  or  //* 
projection.  Initial  conditions  for  p  >0  trivially  follow  from  the  solution  at  the  end  of 
the  previous  time  step. 

A  finite  element  solution  of  Eq.  (2)  is  obtained  by  approximating  by  a  finite 
dimensional  subspace  K  of  piecewise  bilinear  polynomials  on  T.  The  finite  element 
solution  U  satisfies 


(V.U,)  +  (V,f(-,-,r,U,U„U,))  +  A(V.U)  =  J  y'^gsds, 

3rp)3£ljy 

for  all  Ve /(To-  (3a) 

P(Uo),  p  =  0 

The  projection  P  at  p  =  0  is  obtained  by  constructing  a  piecewise  bilinear  approxima¬ 
tion  of  Ufl.  For  p  >  0,  we  proceed  in  a  similar  manner  except  that  we  construct  inter- 
polants  using  the  finest  grid  solution  available  at  r  =  p". 


Temporal  discretization  of  Eq.  (3)  is  performed  by  the  backward  Euler  method; 
thus,  we  determine  U'^U.y)  as  discrete  approximation  of  U0c,y,q)  by  solving 


(V.U«)  +  Ar  [(V.f(-.-,r.U^Ui.U/))  +  A  (V,U^)]  =  (V,UP)  + 

^  I  y'^gN(x^^)ds,  forallVe^0»  W 

dIT 

Initial  conditions  for  the  discrete  system  in  Eq.  (4)  follow  the  lines  of  Eq.  (3b)  for  the 
senu-discrete  system. 

A  posteriori  estimates  of  the  discretization  error  of  the  solution  of  Eq.  (4)  are 
obtained  by  means  of  a  p-refinement  technique.  To  begin,  we  calculate  a  second  solu¬ 
tion  U^(x,y)  of  Eq.  (2)  using  piecewise  quadratic  polynomials  in  space  and  tra¬ 
pezoidal  rule  integration  in  time.  This  solution  is  higher  order  in  both  space  and  time 
than  the  solution  of  Eq.  (4);  thus,  the  difference  ||U'^  -  furnishes  an  estimate  of 
the  discretization  error  of  .  The  computational  efficiency  of  this  procedure  can  be 
substantially  improved  by  using  the  nodal  superconvergence  property  of  finite  element 
methods  for  parabolic  problems  [1,  10].  Nodal  superconvergence  implies  that  bilinear 
finite  element  solutions  converge  at  a  faster  rate  in  space  at  nodes  than  elsewhere. 
These  considerations  imply  that  can  be  calculated  as 

+  EHx,y),  (5) 

where  (x  O’ )  is  a  piecewise  bilinear  function  and  (x  ,y )  is  a  piecewise  serendipity 
function  (a  biquadratic  polynomial  less  a  quartic  term)  that  vanishes  at  the  nodes  of  T . 
Specifically,  we  find  that  U'^(x,y)  satisfies 


(V, 


O’?  -  UP 
At 


)  +  •A[(V,f(•,•,q,C^0,^U;))  +  (V.f(-,-,p,UP,UP,UP))]  + 


V 


'I.*' f ,*  •.•'•-I 
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‘/i[i4(V,U^)  +  i4(V,lIP)]  = '/4  I  V^g^(x,y,^)^&  + '/4  J  gsO^^y  ^P)<^ 

dr  p^dOw  dr  pdf2w 

forallVeATo-  (6) 

Thus,  a  trapezoidal  rule  integration  step  is  performed  using  the  backward  Euler  solu¬ 
tion  UP  (jc  ,y )  as  an  initial  condition.  Both  Eqs.  (4)  and  (6)  are  a  nonlinear  algebraic 
system  which  we  solve  by  Newton’s  method.  In  order  to  reduce  the  computational 
effort  associated  with  assembling  and  solving  Eq.  (6),  the  Jacobian  of  Eq.  (4)  is  used 
for  both  Newton  iterations.  The  solution  of  Eq.  (4)  is  obtained  first  and  the  result 
U^  (x  O' )  is  used  as  an  initial  guess  for  (z  ,y ). 

The  piecewise  quadratic  correction  (x  O' )  satisfies 

(V,[(CJ^ -»-£'?)  -  (UP-t-EP)]/A/)  +  H[(V.f(-,-,q,U^+E^U;j^+E’,U;+E/)) 

+  .UP-i-EP  .UjP-t-EjP.UP-f-Ep)]  -1-  VziA  (V,U'^-hE^)  +  /I  (V,UP-i-EP)] 

= J  V^g^(x,y,q)<is  +  Vi  J  V^g^(x,y  for  all  Ve  (7) 

W"  pdOw  dr  pd£Jw 

As  noted,  the  space  consists  of  piecewise  serendipity  functions  that  vanish  at 
the  vertices  of  the  elements.  Trivial  initial  conditions  arc  used  in  the  solution  of  Eq. 
C7)  for  p  >  0.  Interpolated  values  of  the  initial  error  u°(x  j )  -  Ufx  ,y  ,0)  onto  are 
used  at  p  =0. 

Linear  systems  associated  with  the  application  of  Newton’s  iteration  to  Eqs.  (4), 
(6),  and  (7)  are  solved  by  the  Lanczos  acceleration  of  the  Jacobi  iterative  method  as 
implemented  in  the  iterative  solution  package  FTPACK  of  Young  and  Mai  [12]. 


*. 


IV.  EXAMPLES. 


We  consider  a  sequence  of  three  linear  problems  that  are  designed  to  illustrate  the 
performance  of  our  enor  estimation  and  local  refinement  procedures  and  convergence 
of  the  Schwarz  iteration.  Our  results  are  very  preliminary,  and  additional  computa¬ 
tional  work  and  analysis  will  be  necessary  before  firm  conclusions  can  be  drawn. 

Performance  of  our  error  estimation  technique  is  measured  by  the  effectivity  ratio 

.  (8) 

l|u(x.y.,7)  -  U^lh 

which  is  a  ratio  of  the  estimated  to  the  actual  error  in  the  //^  norm.  Ideally,  the 
effectivity  ratio  should  approach  unity  as  the  mesh  is  refined  and  should  not  differ  sub¬ 
stantially  from  unity  over  a  large  range  of  mesh  spacings.  The  convergence  of  our 
error  estimate  to  the  true  discretization  error  has  been  established  for  one-dimensional 
linear  problems  [10]. 

Example  1.  Consider  the  linear  constant  coefficient  heat  conduction  problem  on 
Q  :=  {  (j:  O' )  I  0  <  X  ,y  <  Jt } 


u,  =  +  Uyy),  (xo')e  fi,  r  >  0,  (9a) 

u(xo',0)  =  sinxsiny,  (x,y)6  (9b) 

u(;x,y,t)  =  0,  (x,y)e9n,  f  >  0.  (9c) 

The  exact  solution  of  this  problem  is 

uCxo'.O  =  e"'«(x,y,0).  (10) 


We  solved  Eq.  (9)  for  a  single  time  step  on  uniform  grids  having  equal  temporal 
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and  spatial  mesh  spacings  of  it//,  /  =  10,  20,  40.  The  exact  enor  and  effectivity  ratio 
are  presented  in  Table  1.  The  results  indicate  that  the  finite  element  solution  is  con¬ 
verging  at  a  linear  rate  and  that  the  effectivity  ratio  is  converging  to  unity. 


7 

||u(x.y.Ar)  - 

6 

10 

0.1578 

1.050 

20 

0.0882 

1.012 

40 

0.0469 

1.003 

Table  1.  Error  and  effectivity  ratio  for  one  time  step  and  uniform  spatial 
meshes  of  spacing  7t/7  for  Example  1. 

Example  2.  Consider  the  forced  heat  conduction  equation  on 
{ (xo')  I  0  <x,y  <  1 } 

Ur +/(xo'.0  =  u„  +  Uyy,  (xo')en,  t  >  to,  (11) 

with  f(x,y,t)  and  the  initial  and  Dirichlet  boundary  conditions  specified  so  that  the 
exact  solution  is 


u{x,y,t)  =  (12) 

With  to  =  0.5,  we  solve  Eq.  (11)  for  one  time  step  on  uniform  grids  having  equal  tem¬ 
poral  and  spatial  meshes  of  1/7,7  =  10,  20,  40.  Results  similar  to  those  of  Example  1 
are  displayed  in  Table  2.  Thus,  once  again,  the  error  is  converging  to  zero  at  a  linear 
rate  and  the  effectivity  ratio  is  tending  to  unity  and  is  close  to  unity  for  all  meshes.  In 
this  example,  as  opposed  to  Example  1,  the  effectivity  ratio  appears  to  be  converging 
to  unity  from  below.  In  practice,  an  upper  bound  is  more  suited  to  an  adaptive  local 
refinement  procedure. 
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J 

l|M(x,v.Ar)-f;^ll, 

6 

10 

0.6796 

0.996 

20 

0.3383 

0.998 

40 

0.1668 

0.999 

Table  2.  Enx>r  and  effectivity  ratio  for  one  time  step  and  uniform  spatial 
meshes  of  spacing  nU  for  Example  2. 

We  also  solve  Eq.  (11)  for  0  =  t  ^  1  using  the  adaptive  local  refinement 
strategy  of  Section  II  with  a  tolerance  of  0.05  and  an  initial  10  x  10  mesh  having  a 
time  step  of  0.1.  Surface  renditions  and  contour  plots  of  the  solution  at  r  =  0.3,  0.5, 
and  0.8  are  shown  in  Figures  4  and  5,  respectively. 

Example  3.  Consider  the  forced  heat  conduction  equation  (11)  on 
:=  { (J^O’)  I  0  <  ac,y  <  1  )  with  f  (x;yj)  and  the  initial  and  Dirichlet  boundary  con¬ 
ditions  specified  so  that  the  exact  solution  is 

u(x^,t)  =  1.0  -  tanh[10(x+y-f-0.45)].  (13) 

This  example  is  used  to  verify  convergence  of  the  Schwarz  alternating  principle.  The 
problem  is  solved  for  a  single  time  step  with  Iq  =  0.5  on  an  initial  uniform  coarse 
10x10  mesh  having  a  time  step  of  0.1  and  a  tolerance  of  0.05.  Refinement  was 
needed  at  the  initial  time  and  10  local  grids,  as  shown  in  Figure  6,  were  introduced. 
The  initial  coarse  mesh  is  also  shown  as  a  reference.  Schwarz  iterations  were  per¬ 
formed  on  these  grids  and  we  measure  the  difference  in  successive  solutions  on  alter¬ 
nating  grids  on  the  portions  of  the  boundaries  of  each  local  grid  in  regions  where  they 
overlap.  The  maximum  such  difference  after  each  Schwarz  iteration  is  shown  in  Table 


^:aagK  sesssae^ 


3.  It  appears  that  the  iteration  is  converging  at  nearly  a  quadratic  rate. 


& 


Iteration  Maximum  Difference 

1  _ 0.1506 

2  _ 0.0114 

3  _ 0.0016 

4  _ 0.0004 


Table  3.  Maximum  difference  between  solutions  on  the  boundaries  of  over¬ 
lapping  grids  after  each  Schwarz  iteration. 


V.  CONCLUSIONS. 


We  developed  an  adaptive  local  mesh  refinement  procedure  for  nonlinear  para¬ 
bolic  systems  on  rectangular  regions.  A  complex  tree  data  structure  is  used  to  manage 
a  nest  of  local  overlapping  grids.  An  implicit  finite  element  solution  strategy  using 
piecewise  linear  approximations  and  the  backward  Euler  method  is  formulated.  We 
obtain  an  estimate  of  the  local  discretization  error  of  these  finite  element  solutions 
using  a  p-hierarchical  approach  with  piecewise  serendipity  approximations  and  tra¬ 
pezoidal  rule  integration.  The  Schwarz  alternating  principle  is  used  to  calculate  boun¬ 
dary  conditions  on  portions  of  local  grids  that  overlap. 

Our  results  indicate  that  the  error  estimation  procedure  converges  to  the  exact 
local  error  as  the  mesh  is  refined.  As  noted,  a  proof  of  this  convergence  has  been 
established  for  certain  linear  one-dimensional  problems  (cf.  Moore  and  Flaherty  [10]). 
It  should  be  possible  to  construct  a  proof  of  convergence  of  the  two-dimensional  eiror 
estimate  using  the  ideas  developed  in  the  one-dimensional  case.  The  use  of  the 
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Schwarz  alternating  principle  also  appears  to  be  a  very  efficient  method  of  calculating 
boundary  conditions  in  overlapping-grid  regions. 

We  arc  encouraged  by  the  performance  of  our  methods  on  these  preliminary 
problems;  however,  several  aspects  of  our  approach  need  improvement.  The  Lanczos 
iteration  used  to  solve  the  linear  system  appeared  to  be  far  less  than  optimal.  The 
stopping  criteria  used  in  the  ITPACK  [12]  implementation  was  too  conservative  for 
our  applications.  Creation  of  local  solution  grids  is  difficult  and  complex  near  domain 
boundaries.  At  present  we  know  of  no  way  of  improving  this  defect.  We  have  plans 
of  extending  our  methods  to  non-rectangular  domains  using  an  overlapping-grid  mesh 
generation  procedure. 
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